arXiv:1507.05182vl [math.NA] 18Jul2015 


An asymptotic preserving scheme for 
kinetic models for chemotaxis phenomena 

A. Bellouquid" and J. Tagoudjeu^ 

“ Cadi Ayyad University, ENSA Marrakech, Morocco 
a.bellonquid@uca.ma 

^ University of YaoundAT I, ENSP and CETIC, Cameroon 
jtagoudjeu@gmail.com, jacques.tagoudjeu@polytechnique.cm 


Abstract 


In this paper, we propose a numerical scheme to solve the kinetic 
model for chemotaxis phenomena. Formally, this scheme is shown to be 
uniformly stable with respect to the small parameter, consistent with the 
fluid-diffusion limit (Keller-Segel model). Our approach is based on the 
micro-macro decomposition which leads to an equivalent formulation of 
the kinetic model that couples a kinetic equation with macroscopic ones. 

This method is validated with various test cases and compared to other 
standard methods. 
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1 Introduction 

Chemotaxis is a process by which cells change their state of movement react¬ 
ing to the presence of a chemical substance, approaching chemically favorable 
environments and avoiding unfavorable ones. In the simple situation where we 
only consider cells and a chemical substance (the chemo-attractant), a model 
for the space and time evolution of the density n = n{t, x) of cells and the chem¬ 
ical concentration S = S{t,x) at time t and position x has been introduced by 
Patlak [U] and Keller-Segel [301 reads: 


^ + V, • (nx(5) ■ V,5 - Dr, ■ V,n) = 0, 
§-DsAS = H{n,S), 


( 1 ) 


where denotes the gradient with respect to the spatial variable, and the 
positive constants Ds and are the diffusivity of the chemo-attractant and 
of the cells, respectively, and x is the chemotactic sensitivity. 
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In general the substance S does not only diffuse in the substrate, but it can 
also be produced by the bacteria themselves. The role of the function H{n, S) 
is to describe the interaction between both quantities. One typical example is 
given by 

H{n, S) = an — bS, (2) 

which describes the production of the chemo-attractant by the bacteria at a 
constant rate a as well as chemical decay with relaxation time ^. Since the bac¬ 
terial movement is directed toward the higher concentrations of S, the coupling 
is attractive. A deep insight into the phenomenological derivation of Keller and 
Segel types models is given in the survey |23] . 

The behavior of this system is now quite wellknown: in the one-dimensional 
case, the solution is always global in time. In several space dimensions, if initial 
data are small enough in some norms, the solution will be global in time and 
rapidly decaying in time; while on the opposite, it will explode in finite time at 
least for some large initial data. The simplicity, the analytical tractability, and 
the capacity to replicate some of the key behaviors of chemotactic populations 
are the main reasons of the success of this model of chemotaxis. In particular, the 
ability to display auto-aggregation, has led to its prominence as a mechanism for 
self-organization of biological systems. Moreover, there exists a lot of variations 
of system o to describe biological processes in which chemotaxis is involved. 
They differ in the functional forms of the three main mechanisms involved in the 
chemotactical movement. They are: the sensing of the chemoattractant, which 
has an effect on the oriented movement of the species, the production of the 
chemoattractant by a mobile species or an external source, and the degradation 
of the chemoattractant by a mobile species or an external effect. The surveys 
|24| and [5] provide a detailed review and critical analysis of the qualitative 
properties of the solutions to problems related to the application of Keller and 
Segel models to various biological contexts. 

Another point of view was introduced by a mesoscopic description of these 
phenomena bridging from stochastic interacting particle systems to macroscopic 
equations. This middle ground consists in describing the movement of cells 
by a "run & tumble" process [33 [33 • The cells move along a straight line 
in the running phase and make reorientation as a reaction to the surrounding 
chemicals during the tumbling phase. This is the typical behavior that has been 
observed in experiments. The resulting nondimensionalized kinetic equation, 
with parabolic rescaling, reads 

' e§+vV,f = ^r(SJ), 

< §-DsAS = H(n,S), (3) 


. f{0,x,v) = fo{x,v), 

where f(t,x,v) denotes the density of cells, depending on time t, position x € 
O C and velocity v € V C T is an operator, which models the change of 
direction of cells and e is a time scale which here refers to the turning frequency. 


2 




The function S{t,x) is the chemical concentration, where n denotes the density 
of cells, and is given by 



(4) 


Starting with the kinetic equation ([3]), one can (at least formally) derive 
the macroscopic limit ([1]) as e —>■ 0. The details will be described in Section 3. 
Other asymptotic limits, such as the hyperbolic limits have been investigated 


in Refs [mil mill [22]. 


The general context of this paper is the development of numerical schemes 
for solving the kinetic equation that are uniformly stable along the transition 
from kinetic regime to the fluid regime. The main difficulty is due to the term 
- which becomes stiff when e is close to zero (macroscopic regime). In this case, 
solving the kinetic equation by a standard explicit numerical scheme requires 
the use of a time step of the order of e, which leads to very expensive numerical 
computations for small e. To avoid this difficulty, it is necessary to use an 
implicit or semi-implicit time discretization for the collision part. In fact, such 
numerical schemes should also have a correct asymptotic behavior: for small 
parameter e, the schemes should degenerate into a good approximation of the 
asymptotics (Keller-Segel model) of the kinetic equation. This property is often 
called "asymptotic preserving", and has been introduced in |26| for numerical 
schemes that are stable with respect to a small parameter e and degenerate into 
a consistent numerical scheme for the limit model when e ^ 0. 

Considering that this paper deals with asymptotic preserving scheme (AP), 
one also has to mention that there are different approaches to construct such 
schemes for kinetic models in various contexts. We mention for instance works 
based on domain decompositions, separating the macroscopic (fluid) domain 
from the microscopic (kinetic) one (see [141 (Hj). There is another kind of 
(AP) schemes for kinetic equations, which are based on the use of time relaxed 
techniques where the Boltzmann collision operator is discretized by a spectral 
or a Monte-Carlo method (see [201 l38l [39l [40]). Other techniques have also 
been developed to design multiscale numerical methods which are based on 
splitting strategy [121 I2Z1 UHl ISTj . penalization procedure dZl [IH1I2S1I13] or 
micro-macro decomposition which first was used by Liu and Yu for theoretical 
study of the fluid limit of the Boltzmann equation [34]. It was then used to 


develop an AP scheme for different asymptotics (diffusion, fluid, high-held, ...), 
see |Z||8l|9l[l3lll5ll25l|321l33]. 


In this work, we extend a method of micro-macro decomposition in order to 
construct asymptotic preserving schemes (AP) for kinetic equations describing 
chemotaxis phenomena. Our strategy consists in rewriting the kinetic equa¬ 
tion as a coupled system of kinetic part and macroscopic one, by using the 
micro-macro decomposition of the distribution function. Indeed, this function 
is decomposed into its corresponding equilibrium distribution plus the devia¬ 
tion. By using a classical projection technique, we obtain an evolution equation 
for the macroscopic parameters of the equilibrium coupled to a kinetic equation 
for the non-equilibrium part. Although our approach is rather general to apply 
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to a very large class of collision operators, the numerical tests shown in our work 
were obtained with very simple model. 

The outline of this paper is the following. In Sect. 2, we present the kinetic 
model and its properties. The micro-macro decomposition, the corresponding 
formulation of the kinetic equation, and the macroscopic limit are presented in 
Sect. 3. Our numerical scheme is presented in Sect. 4. Finally, our method is 
demonstrated with a numerical tests in Sect. 5. 


2 The kinetic model 

In this section we briefly describe the kinetic models we used. The turning 
kernel 'T in the kinetic equation (|31) needs to be specified. The turning kernel, 
7” measures the probability of velocity jump of cells from v' to v. To derive the 
Keller- Segel equation m as £ —0, one has to incorporate both 0(1) and 0(e) 
scale into T. In the following work, as in [D m mi [37] we assume herein that 
the turning operator is of the form 

T{SJ) = To{f) + eTiiS){f), (5) 

where To represents the dominant part of the turning kernel modeling the tum¬ 
ble process in the absence of chemical substance, supposed independent of S 
and 7i(S') is the perturbation due to chemical cues. We first mention some 
assumptions on the turning operators To and Ti(S'): 

• The operators To and 7i preserve the local mass: 

f %{f)dv= [ ri{SJ)dv = 0, (6) 

Jv Jv 


for any S' > 0. 


• There exists a bounded velocity distribution M{v) > 0, independent of x 
and t, such that: 


1. The flow produce by the equilibrium distribution M vanishes, and M 
is normalized: 


/ vM{v)dv = 0, / M{v)d- 

Jv Jv 


'v = 1. 


2. The detailed balance 


(7) 


Toiv\v)M{v) = Toiv,v')M{v') 


( 8 ) 


holds. 

3. The kernel Tq{v, v') is bounded, and there exists a constant cr > 0 such 
that 

Ti^{v,v') > crM, € V X V, X G t > 0. (9) 
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The most commonly used assumption on the turning operators 71 , * = 0,1, is 
that they are both linear integral operators with respect to / and read: 


Tz{S,f)= f {T,{S,v,v')f{t,x,v')-Ti{S,v',v)f{t,x,v))dv'. (10) 

Jv 

Here the turning kernel Ti{S,v,v') describes the reorientation of cells, i.e. the 
random velocity changes from v to v' and may depend on the chemo-attractant 
concentration and its derivatives. 

Technical calculations (see [TJ HT]), namely by integration over v, inter¬ 
changing V hy v'^ and using (l8|), yields the following equality: 


To ( 5 ) 


Hv) 

M{v) 


dv 


5 XX*'"' 


gW) 9{v) 


M{v') 
h{v') 


\M{v) M{v') 


M{v) 
dv dv', 


( 11 ) 


where 

'^[M] = - [To{v,v')M{v') +To{v',v)M{v)). 

In particular Eg. lfTTl) shows that the operators To, is a self-adjoint and the 
following equality: 


Toih) 


Kv) 

M{v) 




hW) 

M{v') 


2 

dvdv' > 0 


( 12 ) 


holds true. 

Moreover, for jy h{v) dv = 0, Ea. dT^ and the estimate ([H]) yield the follow¬ 
ing inequality: 


- I Toih) 
Jv 


M{v) 


dv > 


> 



h{v') 

M,{v') 


2 

dv dv' 


(13) 


which shows that To is a Fredholm operator in the space L'^iV, Therefore, 

the following result defines the properties of the operator To' 


Lemma 1. Suppose that Assumptions (Q)-® hold. Then, the following prop¬ 
erties of the operators To hold true: 

i) The operator To is self-adjoint in the space 



ii) For f G Lf, the equation To{g) 
which satisfies 


f, has a unique solution g G Lfi 



/ g{v) dv = 0 if and only if / f{v) dv = 0. 

Jv Jv 
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in) The equation Toig) = v M(v), has a unique solution that we call 0(v). 
iv) The kernel of % is N{To) = vect{M{v)). 


3 The micro-macro decomposition of the kinetic 
model and macroscopic limit 

3.1 The micro-macro decomposition 

Let {f,S) be a solution of (|31). We decompose / as follows 

/ = M{v)n + eg, 

where n is the density given by (|3]). Then (g) = 0, where {g) = Jy gdv, and one 
has: 

d{Mn) dg 1 ^ ^ , 

—— + e-£ + -vM-\7,cn + v\/xg=-To{g) 

at at e e 

+ \Ti{S){M{v)n)TT,\S){g). (14) 

Now, we use a projection technique to separate the macroscopic and micro¬ 
scopic quantities n(t, x) and g(t, x, v). Moreover, let Pm, denotes the orthogonal 
projection onto iV(7o)- Then 

Pnih) = {h)M, for any h G L‘^{V, 

M (v) 

so that one has the following: 

Lemma 2. One has the following properties for the projection Pm- 

{I-PM){Mn)=PMig) = 0, 

(/ - Pm){vM • V^n) = vM ■ Va,n, 

(/ - PM){Ti{S){M[v)n) = Ti{S){M{v)n), 
{I-PM){ri{S){g)) = Ti{S){g). 

Proof. The hrsts equalities are trivials since Pm{M) = M, and (g) = 0. As the 
flux produced by M is zero then (/ — Pm){vM ■ Va,n) = vM ■ Vxn. 

By using dH]), one deduce that PM{Ti{S){h)) = 0 for any h G L?, and then the 
third and the fourth equality are completed. □ 

Taking the operator I — Pm into the equation (fl^ and using Lemma [2l 
one deduce 

+ -vM ■ V,n + (/ - Pm){v ■ V^g) = Wo{g) + Wi{S){M{v)n) + Ti{S){g). 
ot e £ £ 
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Integrating (ICT) over v, yields 


dn 

'm 


+ {v ■ '^xg) = 0 . 


(16) 


The micro-macro formulation finally reads: 

' e|f + \vM ■ V,n + (/ - PM){y ■ V^g) = \%{g) + \Ti[S){M{v)n) 

+r,{S){g), 

^ + {n-V,g)=0, 

, ^-Ds^S = H{n,S). 


Equations (1171) are the micro-macro formulation of the kinetic equation ([2]) that 
we want to use to design our AP scheme. The following proposition shows that 
this formulation is indeed equivalent to the kinetic equation ([S]) . 


Proposition 1. 

1. Let {f,S) be a solution of (0) with initial data (/o,S'o)- Then (n,g,S), 
where n = (/) and g = -(/ — M{v)n) is a solution to a coupled system 
m with the associated initial data: 


n{t = 0) = no = (/o), 

g{t = 0) = go = ^{fo - M{v)no), (18) 

S{t = 0) = So. 


2. Conversely, if (n, g, S) satisfies system 1171) with initial data (no, go, So) 
such that (go) = 0 then f = M{v)n + eg is a solution to kinetic model (0) 
with initial data fo = M(v)no -\- ego, and we have n = (/), and (g) = 0. 

Proof. The proof of i) is detailed above. For ii) consider (n, g, S) solution of 
di). We set / = M{v)n + eg and we show that / is a solution of kinetic model 
©■ From dD, one has 

f) f 1 11 

- M{v)— + -vM .y,n + v - Pm{v ■ V^g) = ^To{f) + -ri{S){f), 

and then 

- M{v)^ + -V ■ V,/ - M{v){v ■ V.ff) = 4ro(/) + Wi{S){f). 
ot ot e e 

Therefore using dl), one obtains ©. The property (g) = 0 is obtained by 
integrating (fTK)) over v, using © and the property of the initial data. This 
completes the proof. □ 
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3.2 The macroscopic limit 

In this subsection, the formal derivation of the macroscopic model is performed 
starting from the meso-macro model ([3]). The macroscopic model has been de¬ 
rived mathematically in mi. Hereafter, we will see that the formal derivation 
is really straightforward starting from (fT71) (compared to the equivalent for¬ 
mulation of ([31)), since the micro-macro model is well suited to deal with the 
asymptotic model in the diffusion limit. Indeed for small e, the first equation 
of (fT71) by using ([HI) and ([7]) gives 

g = %-\vM ■ V,n) - %-\ri{S){M{v)n)) + 0{e). (19) 


Inserting (IT31) into yields the asymptotic model (coupled with the concen¬ 
tration equation for S'): 


t - {v.V,{%-Hri{S){M{v)n))) = 0{e), 

§-DsAS = H{n,S). 


( 20 ) 


Using hi) of Lemma [U one has as To is self adjoint operator in 
the following: 


{v.S/,{%-\T,{S){M{v)n))) 


(roiOiv)). ^ (V^Ti (S) (M(^)n)) ^ 
divx (^-^^nri{S){M{v))^ , 


and then the macroscopic model (1201) becomes 


^ -I- div^ (na(S) - ■ V^n) = 0{e), 

§-DsAS = H{n,S), 


( 21 ) 


where Dn and a{S) are given by 


Dn = - [ V(g>9{v)dv, aiS) = - f ^P-ri{S){M){v)dv. (22) 

Jv Jv M{v) 

The approach we have developed is quite general. The Keller-Segel model 
is now given. 


3.3 The Keller-Segel model 

Let us first consider the following task for the probability kernels: 

To{v,v') = oM{v), cr > 0. 

Consequently, the leading turning operators To become relaxation operators: 

%{g) = -cj(g-{g)M). (23) 
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In particular, 6 and the diffusion tensor are given by: 

1 If 

0 {v) = - M{v), Dn = — / v®vM{v)dv. 

O' (T Jv 

Moreover a{S) is given by: 

a{S) = - [ vTi[S]{M{v))dv. 

O' Jv 

The relaxation kernels presented in together with the choice 

T,[S]=Ks{v,v’)-V,S, 

where Ks{v,v') is a vector valued function, leads to the model 

ri[s]{M) = h{v,s)-v^s, 

where 

h{v,S) = j (^Ks{v,v')M{v') — Ks{v',v)M{v)^dv'. 
Finally, the function a{S) in (l25ll is given by 

a(5) = x{S) ■ 

where the chemotactic sensitivity x(S') is given by the matrix 

x(S') = — / V 1^1 h{v,S)dv. 

O' Jv 


( 24 ) 


(25) 


(26) 


Therefore, the drift term divx{n Q:{S)) that appears in the macroscopic case 
stated by (EOI) becomes: 

diVx{na{S)) = div^ (nxiS) ■ V^S ), 

which gives a Keller-Segel type model m- 


^ + divx (n x[S) -VxS - Dn-Vxu) = 0{e), 


f-DsAS = H{n,S). 


(27) 


4 Numerical methods 

Here in this section, we consider Problem o, subject to the following initial 
conditions: 

f{0,x,v) = fo{x,v) and S{0,x) = So{x). (28) 
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It has been shown that problem ([3]) is equivalent to the following micro-macro 
formulation: 


^ + -^vM ■ VxU. -|- i(/ — Pm){v ■ Vxff) — -^To{g) 

+^TiiS)iM{v)n) + lTi{S){g), 

< (29) 

^ + {v V.g) = 0, 

^ §-DsAS = H{n,S), 

subject to the following initial conditions: 

n{t = 0) = no = (/o), g{t = 0) = go = -(/o - M{v)no), S{0,x) = S'o(x). (30) 

e 

The discretization of problem (I29I) - (I301) is carried out with respect to each 
independent variable (time, velocity and space). 

4.1 Time discretization 


The treatment of the time variable of problem (I331) - (I301) can be done by using 
varieties of methods such as finite difference and variational methods. Finite- 
differencing the derivative in time is the widely used approach. 

The time interval [0,r] is divided into N times steps as follows: 


to = 0, tk+i = tfc -I- At, 0 < fc < iV, 


where At = -^ is the time step. The approximation of n{t,x) and g{t,x,v) at 
the time step t^ are denoted respectively by « n{tk,x) and g^ ~ g{tk,x,v). 
Using an implicite scheme for the stiff term -irToig) ^md an explicite for the 
other terms in the first equation in (I29L one obtains : 




At 


= -^uM-V,n^'--(/-PM)(«-V,/) 

e 


+ -TiiS'^Kg^). 

e 

Substituting g by in the second equation of ([33)) yields 

ijq ^“1“ 1 nrj k 

--+ {v ■ V.g'^+i) = 0. 

Replacing n in the third equation by one has: 

^fe+1 _ i^k 


At 


-DskAS'‘+^ = il(n'=+\5''=+i). 


(31) 


(32) 


(33) 


10 







Proposition 2. The time discretization of the first and second equa¬ 
tion of system is consistent with the first equation of system i2CH) when 

e —^ 0. 

Proof Formally, we have from (IHTl) : 

5 "+' = ff" + ^ [Ti (5'=) {M{v)n>^ + eg'^) - vM ■ 

-e{I-PM)v-V:,g'^]. (34) 

Since the operator —To is self adjoint and positive, (/ — -^To) is self-adjoint 
and positive definite thus invertible for At > 0. Therefore one has 

[g^ + ^[Ti{S^){M{v)n^+eg^)-vM-V,n'^ 
-e{I-PM)v-V^g^]). (35) 

Developing the right hand side of (l35)l with regard to e when e —>■ 0, yields: 

/+i = [vM • V,n^] - [Ti{S^){M{v)n^)\ + 0{e). (36) 

Substituting g^^^ in (15^ leads to 

I ]_ 

"" aP [uM-V, n'=])-(u-V,ro-' [TiiS'^){M{v)n'^)]) = 0{e), 

(37) 

which is consistent with the first equation of system (EUl) when e —?> 0. □ 


4.2 Spatial discretization and velocity dicretization 

In the following, we present the methods in the case of 1-dimensional spatial and 
velocity discretization. The phase-space interval is denoted by [XmimXmax] x 

[Xmiri: Xtjiaxl^ where Vmin — Xmax ^ 0 - 


4.2.1 Spatial discretization 

For the spatial discretization, a difference method based on control volume ap¬ 
proach and cell averaging is used. The numerical grid is defined by: 


R^x — { 2:2 ,^, 0 ^ ^ ^ 


^max ^min 

Ax 




(38) 


where Ax > 0 is the spatial mesh size, Xq = Xmin, Xi = Xi-i + Ax {1 < i < N^) 
and x^_^_l = (xi+i-|-Xi)/2 (0 < i < —1) are the cell center points. Proceeding 

as in [H], the microscopic equation (1311) is discretized at points Xj_|_i while the 
macroscopic equation (1321) and the diffusion equation (1331) are discretized at 
points Xi- The approximation of n(t,x), g{t,x^v) and S{t,x) at the considered 
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spatial points and at the time step tk are denoted by k, n{tk,Xi), i (v) ~ 

2+ 2 

g{tk,x^_^_l,v) and « S{tk,Xi) respectively. 

Setting v~^ = max(0,r’), v~ = min(0,r!), the following spatially discrete 
forms are obtained: 


3i+i, ^i+h 1 




At 

1 I 


+ -(/-Pm) 


Ax 


+ V 




_ ^i+ 


Ax 


_ 71^ 


Ax 




At 


k+l k+1 , 


9^+ 


Ax 


= 0 , 


(39) 

(40) 


and 


qk+l Qk 

7 -) 

-1- jyqk - 

At 


^fc+1 ~h 


''i+i 


(Ax)2 


= Hinr\Sr^). 


(41) 


Proposition 3. From the discretization (E3)-|2F °f yields the following 
numerical scheme when e —>■ 0.' 


- nf 
At 


+ 


1 


Ax 

- % 

1 

Ax 


-1 


vM ■ 


n^+l - rii 


Ax 


-1 I yM ■ 

Ax 






= 0 


(42) 


which is consistent with the first equation of system i2(A) . Moreover, the approx¬ 
imation of the diffusion term is second order accurate in space. 

Proof. The quantity g^ i is derived from (15^1) as: 

*+ 9 


„fc+l _ j T 

9^+2. - - —>0 


9\+\ 


-^{1-Pm)\v 


+ ^+2 ^ 2 I ®+ 


Ax 


1 ' 


Ax 




-vM ■ 

Ax 
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It follows that: 


fc+i 


= % 


-1 


vM • 




Ax 




0{e) (43) 


as e —> 0. Replacing in (l40ll by its expression of (l43l) and passing to the 

2 

limit yields the relation Proceeding as in the continuous case, the spatially 
discrete form (1421) is consistent with the first equation of system (Eni). □ 


4.2.2 Velocity discretization 

We consider a uniform velocity grid defined as: 

Vav = {Vj = Vmin + jAu, 0 <J < 7V„}, (44) 


where Av = is the velocity step and W € N* is and odd number. 

The approximation of g{t,x,v) at the spatial points and velocity Vj at 

time step tk is denoted by ^ ~ ^(^fci ^i+i iThe phase-space discrete 
formulation of system dlB writes: 


fc+i 






^ Ax 


ToAgT^) + TiAsLAiMAAlA -• 


q-i-i 


n+i 


Ax 




- A 


At 


+ ( Vi 


k+i 




Ax 


= 0 , 


(45) 

(46) 


and 


Ofc+l Qk 

At 


-D 


sf 


Q^+l eynkAl I cAl + l 

~ ^ = GinA\SAA (47) 


0<k<N,0<i< N^, and 0 < j < iV„. 

In the discrete case, the bracket (.), the projection I - Pmj and the integral 
operators Toj and Tij{SAA defined by: 


Ny-1 

(Gj) = 

Nv-1 

{I - PmAGj) = Gj - M,Av Y, Gk, 

k^O 


(48) 


(49) 
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(50) 


^Nv-1 


Nv-l 


ToAG) = aW ^ To{v,,vi)Gi - G, Y. T^{vuv,) , 


1=0 


1=0 


and 


C N„-1 N^-1 \ 

E T,{SY.,v,,vi)Gi-G, E Ti{SYY^,v,)\ , 

1=0 ^ 1=0 ^ J 

(51) 

where uj = {uj + Uj+i)/2. 

4.2.3 Boundary condition 

For the numerical solution of the kinetic equation ([3]), the following inflow 
boundary conditions are usually prescribe for the distribution function /: 

f{t,Xmin,v) = fl{v), V>Q and f{t,Xmax,v) = fr{v), V < 0, (52) 

which can be rewritten in the micro-macro formulation: 

i{t,xo)Mj + I ^g{t,xi,Vj)+g(t,x_i,Vj)^ = fi{vj), Vj > 0, 

i{t,XNjMj + I (^g{t,XN^+i,Vj) + g{t,Xf^^_i,Vjfj = fr{vj), Vj < 0. 


n( 


The following artificial Neumann boundary conditions are imposed for the other 
velocities m- 


git,x_i,Vj) = g{t,xi,Vj), Vj < 0, 

g{t,XN^+i,Vj) = g{t,xj^^_i,vj), vj > 0 . 

Therefore, the "ghost" points can be computed as follows: 

fc+i f u, >0, 

9-1 i 


fe+i 
^Afx + 


fc-tl 

ds 

A 

o 


f ifrAj) 



J-yx 2 

Vj > 0 . 



(54) 

(55) 

(56) 


It then follows from (l46)l that: 

1+(d*t)) '■;«="S - D ( i«i +< - ”7)4? - 


(57) 


and 






At / 2v~, 


= n 


N, 


■friVj) 


Ax \ e 

-iv,-vi + v-)gYl, 


(58) 
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Additionally, the homogeneous Neumann boundary conditions are pre¬ 
scribed for the concentration S: 


q/c+I C^ + l _ C^ + 1 

0_1 — dlia 

4.2.4 Implementation of the method 


(59) 


To implement our method, we consider H given by ([2). At the time step tk, let 
us define the following vectors: 


M = 


/ A/(uo) \ 


/ Gi \ 

gI 

2 

M{vi) 

, G''' = 

^ M{vnJ ) 


1 ) 


with 




TV'' = 

in 


i ^ 


(St \ 


( F( \ 

n\ 

, = 

Si 

and E'' = 

pk 

K } 


\ 51 ) 


\ Fk / 


with = a{At)n^~^^. 
Setting 


0 AtAv _ 0 AtAv^^^ 

w)i = + To{vi,Vj), 


Nv—1 


1^0 


AtAv 


1=0 

„„ifc AtAv 11 ^ 

= ——Pi."-’ o-"- "■ = 


1,3 


1,3' 


2k _ qlfc fc ... _ 

^1,3,1 J.2 Pi,3^i+i’ ^1,3 


2k _ k I (Af)u; 

2 V 2 1 


Aa; 


then Eq. (HSl) can be recast as: 

(a°-b°)g'"+1 = (bI’^-p++p~+i)g’^^i+b^’^m+p+g':_i-p~g':^3 ( eo ) 

2+ ^ 2-|- 2 ‘2 2 

for 0 < J < Alj, — 1, where 

i^j,i = (1 < J,^ < A„-b 1), 




0 

j-i.o 


(1 < J < A, + 1), 


= 2 

Bh = (1 < J < A„ + 1,2 < Z < AT,), 


dO 


= (1<J< A, + l), 
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Ik 


B 


Ik 


B 




,,lfe ;; Ifc 


, (1 < j < iV„ + 1), 


= <5-1,/-I - {1<J<N, + 1,2<1< N,), 


Ak 


„lfe 


-k X ^Ik 


1,3 


(1 < j <iV,+ 1), 


B 


r)2fc 

r?2k 


2k 

hj,Nv-\-l 


Ak _ A rr'^k 

^i,j-l,0 


{l<j<N^ + 1), 


= <5-1,,-i - <5i/<5_i, (1 < J < iV. + 1, 2 < / < N,), 


w: 


Ik X ^2/e 

,j-l,Ny ~ ^jlV„+lcri,Ar„ 


, (l<j<iV„ + l), 


p± — 
^3 A ~ 

p± — 


P: 


j,N^ + l 


- Mj i{Av)vq ) 


, (l<j<iV„ + l), 


(<5,,, - M,_i(Ai;)< J, (1 < j < + 1, 2 < / < iV,), 


At 


eAx 

^i^3,N.+i - Mj-i{Av)v^J 


(1 < j < A„ + 1). 


Similarly, the vectors and are respectively computed from (I46II 

and (H71) as follows: 


and 

where 


with 


and 


« = <n.i - - Gf_ J, t = 1,2, •. • , A, - 1, 




/ Ofco bko 

bki flfci &fei 


Afc = 


V 


bkN^-1 afcAfx-i bkN^-i 

bkN^ akN^r / 


6fco = — 2co, 

bki — Ci, t — 1, 2, • * • , Aa; 1, 

= ~‘2CNa:, 

Ofci = (1 + 2ci + 6(At)), t = 0,1, • • ■ , Aa; 

{At)Dgk 


Ci = 


{AxyP ’ * — 0, 1, • • • , Aa; 


At Aw /Wo Wat, 




(61) 

(62) 
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The following algorithm is hnally obtained for the numerical solution of the 
Micro-Macro system (ITTl) : Given G°, i, G° i, 5°, TV? , rig, n% . 

For fc = 

1. Solve G^_^i, (j = 0, 1, • • • ,iVa: — 1) from (1^ : 

2+ 2 

2. Compute using (ICTl) : 

3. Compute Ug^^ and using (1571) and (1551) : 

4. Compute G^V and using (|551) : 

2 iVa:+ 2 

5. Solve 5^+^ from (157|) . 


4.3 A time implicit discretization 

The previous dicretization is explicit of the macro part n. It then imposes 
the diffusion restriction on the time step At = 0{{Ax)‘^). To overcome this 
restriction, a time implicit discretization can be applied for the macro part such 
that at the limit, the diffusion term is treated implicitly. Following the idea in 
|15| . a time implicit scheme can be derived the macro part of the micro-macro 
system. It consists to substitute in (1401) by defined as follows: 


~k+i 


where 


= 


-1 ( 
Oe 




n, 


Oe 


I V 


k _ k 

, y,_i 




9r+i 


Ax Ax 


'Hoe = { I -y7o 


and = 

l-\-o 


+ l _ 77 

Ax 


fc+i 


(63) 


The following implicit time discretization of the macro part is obtained: 


At 


^UPg-; I VM- 


Ax 


9 ^ 1 - 9 ^- 1 ' 

’ Ax 


= 0 , 


where 




(64) 
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It can be seen that for small e, 


(vM ■ = -Tq ^ (vM ■ + 0(e) 

and 

3^4= V (ri(S^)(M(v)n^^^)) . 

Trough substitution and using the properties of To, we obtain the following 
discret form of the macro part as e —> 0 


- nf 


At 


-A 


- 2n: 


fe+i 


,fe+i 

b+i 


(Ax)2 




Ax 


= 0 , 


(65) 


which in consistent with a discret form discretization of the macroscopic limit, 
using implicit dicretization of the diffusion term. We remark that there is an 
additional computation of T-Lq^ for the calculation of In the particular case 

where 7o(/) = —<j(f — (f)M(v)), we have from the micro-macro decomposition 
M9i+l) = Therefore 


At 


At / 


+ a At \ 


Ax 


Qk+l _ Afc + 1 

9,+ i 9,_i 


Ax 


= 0 , 


and is obtained by solving the linear system 
(An'^+iAt 


k+i _ ^k+i . 


Ax / ’ 

where A is the tridiagonal matrix A = Tridiag(—a, 2 -|- a, —a) with 

^.k+i _ fc+i At 


a = (AM) 


(At)2(e2 -I- At) 


and + 


( 66 ) 


(67) 


a At 




The boundary conditions are incorporated using (1571) and (1581) . 


5 Numerical results 

We present here some numerical experiments to validate our approach. In our 
tests, the space domain is the interval X = [—1;1], the velocity domain is 
V = [—1; 1]. For all the numerical tests carried out below, the velocity space is 
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divided into iV„ = 64, which can provide good enough accuracy for numerical 
simulations |10| . The equilibrium distribution M{v) and the kernels Tq(v,v') 
and Ti{S,v,v') are set to: 

M{v) = i To{v,v') = M{v), T^{S,v,v') = {v.S/S)+. 

It is clear that Tq and Ti satisfy assumptions (0)-®. The boundary conditions 
are given by f{t,—l,v) = 0,Vu < 0-,f{t,l,v) = 0,Vt! > 0. For the chemoat¬ 
tractant equation, we consider H{n,S) = —S + n and the initial condition 
S(0^x) = 0. The initial condition for the cell density is a gaussian-like peak 
given as: 


n{0,x) = C'Mexp(—80a;^), 

where Cm is a constant determined such that the total mass is M. The initial 
cell distribution function is given as: 

/(0,x,.)=(„(0,x) + =2^)m(,), 

We compare our scheme with: 

• an Explicit-Euler scheme applied to the kinetic equation in the kinetic 
regime and to the diffusion equation in the Keller-Segel limit. 

• an asymptotic preserving scheme obtained from a time splitting method 
applied to the Odd-Even decomposition of the kinetic equation in both 
kinetic and macroscopic regime. 

• Explicit Euler scheme for the kinetic equation 

The simplest time explicit scheme for the kinetic equation reads: 

= /" - ^ K(/" - /ti) + - fh) + (68) 

Because of the stiff term £~^T{S^, this scheme is not asymptotic preserving. 

• Explicit Euler scheme for the Keller-Segel system 

Under the above setting, the Keller-Segel system (EID reads: 

f # = d^{Dd^n - xnd^S), 

\ § = DsAS-S + n, 

associated with homogenous Dirichlet boundary conditions for n (obtained from 
((551) when £ ^ 0) and Neumann boundary conditions for S, where the coeffi¬ 
cients D and x obtained from (l24l) and (l26l) respectively as Dn = ^ and 
X=\- 
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The numerical solution to the first equation of the Keller-Seigel system can be 
computed by using the following standard explicit finite difference method |42| : 


^fe+i 


At 


= D 


H+l 


— 2n^ + n: 


i-l 


(Ax)2 


X- 


(5; 


(c) 




Ax 


(69) 


where 

{l<^<Nx- 1 ), = 0 . 

The diffusion term can also be treated implicitly in (IMll . 

• The Odd-Even parity decomposition asymptotic preserving scheme 


The Odd-Even decomposition asymptotic preserving scheme is based on the 
following parity equivalent form of the initial kinetic equation [iniET]: 


^ + vd^j = ^[{1 + e\vdj;S\)n - {2 + e\dj;S\)r], 
% + = ^[{vdxS)n -{2 + e\dxS\)j], 


(70) 


where the even part r{t,x,v) and the odd part j{t,x,v) of the flux f{t,x,v) 
verify: 

r{t, X, v) = i(/(t, X, v) -k fit, X, -v)), 

(71) 

j{t,x,v) = ^{f{t,x,v) - f{t,x,-v)). 

It then follows that 


fit,x,v) 


r(t, X, v) + sjit, X, —v)), w > 0, 
r(t, X, -v) - ej{t, X, -v)), v < 0, 


(72) 


and 


n = 2 rdv. 

Jo 

Therefore, the vacuum boundary conditions are defined for n > 0 as: 

{r{t,x,v) + ej{t,x,v))\j;=_i = 0, 

{r{t,x,v) - ej{t,x,v))\x=i = 0. 

When e < 1, the parity system cni) can be rewritten as: 

^ + vdxj = ^[{1 + e\vda;S\)n - {2 + e\dxS\)r], 

|l + vd^r = ^[{vdxS)n - {2 + e\dxS\)j] -k (1 - ^)vdxr. 


(73) 


(74) 
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Taking the limit as £ —>■ 0 in (ITU) , yields 

j{t,x,v) = ^{vdoi:S)n - vd^r. 


Replacing j in (1731) and taking into account the homogenous Neumann boundary 
condition prescribed on S, one obtains the following boundary condition: 

( (r - eu5a:r)|x=-i = 0, 


{r + £vda,r%=i = 0. 


(75) 


The odd-even parity asymptotic preserving scheme for the initial kinetic equa¬ 
tion consists of a two steps operator splitting method which combine an implicit 
scheme for the stiff collision term 


dt 


= + £\vdxS\)n - {2 + e\dj;S\)r], 


m = -{2 + £\d^S\)j] -k (1 - ^)vd^r, 


(76) 


at ~ 2e2 

with an explicit scheme for the non stiff transport term 


% + vd,3 = 0 , 

-H vd^r = 0. 


(77) 


For the spatial discretization, a center finite difference can be used for the col¬ 
lision step (1761) and a second order upwind scheme for the transport step (iTTll . 


Numerical tests: In the following, we denote by: 

• MM: the scheme obtained from the micro-macro decomposition, 

• K-S: the scheme for the keller-Segel system, 

• Explicit: the explicit scheme for the kinetic equation, 

• Odd-Even: the odd-even parity asymptotic preserving scheme. 

We have observed that the use of the time implicit discretization for the 
micro-macro model and the Keller-Segel limit give rise to numerical results 
which are very close to those produced by the explicit discretization. Hence, for 
the numerical results presented, we use the implicit approach for the MM and 
K-S schemes. The first test concerns the convergence order of the MM scheme 
computed at time t using the P as : 

/ _ II/Ax(0 - f2Ax{t)\\2 

II/2Ax(0)||2 

where /ax denotes the approximation of / using the spatial grid size Aa;. The 
time step is set to At = . Figure [T] presents the convergence rates obtained 
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with Nx = 80, 160, 320, 640 at time t = 0.1 for e G {1, 0.01, 10“"^, 10“®} and 
the total mass: M = 2tt (Figure [T]). It can be seen that the MM scheme 
converges uniformly since time step does not depend on e. A second order 
convergence is observed in the diffusive regime (e < 10“"^). 



Logio(Ax) 


Figure 1: Convergence order of the method for e G {1, 0.01, 10 10 ®} at 

time t = 0.1 (M = 27r). 


In the following, we set Nx = 200, M = 27r. The time step is set to At = 

at the kinetic regime and At = 0(Aa;) at the diffusive regime (At = Ax/2 
for MM and K-S schemes and At = Ax/40 for the Odd-Even method). We 
mentioned that the comparison of the MM, Odd-Even and Explicit methods, 
the initial condition has been projected at the equilibrium for the last two 
methods at the kinetic regime. We illustrate in Figure [2] the behaviour of the 
MM scheme at different regimes. For different values of e {sk = 2“^, k > 0), 
we plot at time t = 0.5 the density of cells. We also add the result obtained 
with the K-S scheme. It can be seen that the MM scheme is stable as e —>■ 0 
and converges to the Keller-Segel limit. Indeed, for e < 2“^, the profiles of the 
density given by the two schemes are quite the same. 

To check the behaviour of MM scheme in kinetic regime, we compare in 
Figure [3] the density of cells obtained for e = 1 with the MM, Explicit and 
Odd-Even schemes at time t = 0.5. As expected, for both schemes, the results 
are very closed. 

We also illustrate the behaviour of the methods in macroscopic regime. We 
compare in Figure |4] the density of cells obtained for e = 10“® with the MM, 
Odd-Even K-S schemes at time t = 0.5. As expected, for both schemes, the 
results are quite the same. 

We investigate the time evolution of the density using the MM scheme in 
different regimes (e = 1,0.01,10“®). The results are shown in Figures [3 [6] and 
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Figure 2: Density of cells at time t = 0.5 using MM and K-S schemes for e = 2 ^ 
{k G {0,3, 5, 7,9}) and M = 2tt 



Figure 3: Density of cells at time t = 0.5 obtained with MM, Odd-Even and 
Explicit schemes for e = 1 and M = 27r. 
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Figure 4: Density of cells at time t = 0.5 obtained with MM, Odd-Even and 
K-S schemes for e = 10“® and M = 27r. 
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1 


Figure 5: Evolution of the cell Density using MM scheme for £ = 1 and M = 27r. 
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Figure 6: Evolution of the cell Density using MM scheme for e = 0.01 and 
M = 2tt. 



time t 


Position X 


Figure 7: Evolution of the cell Density using MM scheme for e = 10 ® and 
M = 2tt. 
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[71 In each figure, the density seems to evolve to a stationary solution. 

6 Closure looking ahead at research perspectives 

This paper has developed a computational approach to a class of pattern for¬ 
mation models derived from the celebrated Keller-Segel model obtained by the 
underlying description delivered by generalized kinetic theory methods. The 
derivation is based on a decomposition with two scales, namely the microscopic 
and the macroscopic one technically related, as we have seen, by suitable small 
parameters accounting for the time and space dynamics. 

The novelty of our paper is that the computational scheme follows precisely 
the derivation hallmarks by using the same decomposition and parameters. This 
idea improves the stability properties of the solutions with respect to classical 
approaches known in the literature. However, without repeating concepts al¬ 
ready mentioned in the previous sections, we wish to stress that this method can 
contribute to future developments also related to applications. In fact, the need 
of new models in biology is well presented in [53] and |S| to account for a broad 
variety of biological phenomena. Moreover, it is shown in |5] that the so-called 
micro-macro decomposition can lead to an interesting variety of models such as 
models of angiogenesis phenomena. 

Therefore, modeling and computational methods can march together thus 
contribution to a deeper understanding of the specific features of the two dif¬ 
ferent, however, related fields. Indeed, we have in mind not only applications 
in biology, but also to the dynamics of self-propelled particles such as those of 
vehicular traffic as it has been recently shown |5| how macroscopic models can 
be derived from the kinetic description using precisely the micro-macro decom¬ 
position treated in this present paper. 
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